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ABSTRACT 

In  scientific  confutation,  there  is  often  need  for  the  derivatives  as  well  as 
the  values  of  functions  defined  by  computer  programs.  Here,  it  is  shown  how  auto¬ 
matic  differentiation  can  be  carried  out  in  a  modem  computer  language  which  permits 
user-defined  operators  and  data  types.  The  specific  language  used  is  PASCAL-SC,  and 
differentiation  is  implemented  for  variables  of  type  GRADIENT,  which  consists  of  the 
value  of  a  function  of  n  real  variables  and  its  gradient  vector  of  first  partial 
derivatives  with  respect  to  the  independent  variables.  Calculations  of  the  results 
of  operators  or  functions  applied  to  GRADIENT  variables  are  carried  out  according 
to  the  well -known  rules  for  evaluation  and  differentiation  of  sums,  differences,  pro¬ 
ducts,  and  so  on.  Since  the  differentiation  is  performed  at  compile  time,  the  code 
produced  is  comparable  in  compactness  and  execution  time  to  that  obtained  if  numerical 
approximations  are  used  for  derivatives,  and  the  theoretical  and  practical  problems 
associated  with  numerical  differentiation  are  avoided.  PASCAL-SC  source  code  is  given 
for  the  necessary  operators  and  standard  functions,  and  it  is  shown  how  to  prepare 
code  for  arbitrary  differentiable  functions  to  add  to  the  library  if  desired.  The  ef¬ 
fectiveness  of  the  use  of  type  GRADIENT  is  shown  by  an  example  of  the  solution  of  a 
system  of  nonlinear  equations  by  Newton's  method. 
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SIGNIFICANCE  AND  EXPLANATION 


In  scientific  computation,  there  is  often  need  for  the  derivatives  as  well  as 

the  values  of  functions  defined  by  computer  programs.  Examples  include  calculations 

of  sensitivity  and  error  analysis,  unconstrained  and  constrained  optimization,  and 

the  solutions  of  nonlinear  systems  of  equations.  Modem  computer  languages  permit 

user-defined  operators  and  data  types,  so  the  well-known  rules  for  differentiation 

can  be  automated .  Here,  this  is  done  in  the  language  PASCAL-SC  (PASCAL  for  Scientific 

Computation) ,  which  also  offers  the  advantages  of  extremely  accurate  floating-point 

arithmetic  operations  for  real  and  complex  numbers  and  intervals,  and  vectors  and 

matrices  over  these  numerical  data  types.  For  the  purpose  of  differentiation,  a 

data  type  called  GRADIENT  is  introduced.  A  datum  of  type  GRADIENT  is  a  record  of 

the  form  (f  (x)  ,  (f ^  (x)  ,. . .  ,f  (x) ) )  ,  i.e.,  consisting  of  the  value  at  same  point  x  » 

(x, , . . . ,x  )  of  some  function  f  of  n  variables ,  together  with  the  values  at  x  of  its 
l  n 

first  n  partial  derivatives  f^  =  8f/8x^,  i  =  l,...,n.  This  requires  specification 
of  all  allowable  operations  involving  this  data  type.  For  example,  the  multiplication 
operator  *  is  extended  to  the  variables  f  and  g  of  type  GRADIENT  by  the  prescription 
f*g  -  (f (x)*g(x) ,(f(x)*g1(x)+f1(x)*g(x) ,...,f(x)*gn(x)+fn(x)*g(x))) .  That  is  about 
all  there  is  to  it:  The  rules  for  evaluation  and  differentiation  are  applied  to  ex¬ 
pressions  involving  GRADIENT  variables  to  obtain  the  correct  value  and  gradient  vec¬ 
tor  for  the  result.  Since  the  differentiation  is  done  at  compile  time,  the  code  pro¬ 
duced  is  compact  and  efficient,  and  the  practical  and  theoretical  difficulties  as¬ 
sociated  with  numerical  differentiation  are  avoided.  As  shown  by  an  example,  numer¬ 
ical  differentiation  is  inaccurate  and  unstable  even  in  favorable  cases,  and  its  use 
can  render  the  results  of  a  program  meaningless  without  the  user  being  aware  of  the 
loss  of  significance.  Thus,  when  software  for  analytic  differentiation  is  available, 
it  is  foolish  and  perhaps  even  dangerous  to  use  numerical  approximations  to  deriva¬ 
tives. 

PASCAL-SC  source  code  is  given  for  the  necessary  operators  and  basic  standard 
functions.  It  is  shown  how  other  functions  can  be  added  easily  to  the  library,  if 
one  knows  rules  for  their  differentiation  or  how  to  express  them  in  terms  of  differ¬ 
entiable  functions.  An  example  of  the  effectiveness  of  the  use  of  type  GRADIENT  is 
its  use  in  the  solution  of  nonlinear  systems  of  equations  by  Newton's  method,  which 
is  illustrated  by  a  specific  calculation. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive  summary 
lies  with  MRC,  and  not  with  the  author  of  this  report. 


DIFFERENTIATION  IN  PASCAL-SCj  TYPE  GRADIENT 

L.  B.  RAll 

1.  Automatic  differentiation.  In  aciantific  computation ,  there  la  often 
need  for  the  derivatives  as  well  aa  the  valuea  of  the  functiona  being  computed.  For 
problems  of  even  moderate  sire,  however,  it  ia  usually  not  feasible  to  differentiate 
all  expressions  appearing  in  the  program  by  hand;  this  time-consuming  process  can 
also  result  in  additional  errors  which  have  to  be  tracked  down  and  eliminated  from 
the  final  code.  The  alternatives  are  then  to  resort  to  the  use  of  inaccurate  numer¬ 
ical  differentiation,  based  on  difference  quotients,  or  to  have  software  which  will 
enable  the  evaluation  of  the  required  derivatives  automatically  in  the  course  of 
the  computation.  The  latter  approach  is  possible  because  differentiation  proceeds 
according  to  fixed  rules,  and  thus  can  be  automated  [7  j,  fl.2],  [13 1  •  The  imple¬ 

mentation  of  this  method  in  a  modern  scientific  computing  language  is  presented 
here;  a  brief  comparison  of  the  results  of  analytic  with  numerical  differentiation 
will  also  be  given. 

In  order  to  make  an  efficient  implementation  of  the  method  of  automatic  analy¬ 
tic  differentiation,  a  language  is  required  in  which  the  user  can  defined  appropriate 
daJta  typd,  aa  in  ordinary  PASCAL,  and  corresponding  optHttfOKt  on  these  types,  which 
is  a  feature  of  ALGOL-68,  far  example  [15].  The  scientific  computing  language  known 
as  PASCAL-SC  [1],  [19]  has  both  of  these  features,  and  in  addition,  highly  accurate 
arithmetic  based  on  a  general  theory  [  9  1  for  real  and  complex  niasbers  and  intervals , 
as  well  aa  vectors  and  matrices  over  these  numerical  types.  The  latter  capabilities 
are  essential  for  accurate  scientific  computation,  but  do  not  anter  explicitly  into 
the  following  discussion.  Hence,  the  techniques  discussed  below  are  not  limited  to 
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PASCAL-SC,  but  can  be  adapted  readily  to  any  language  in  which  the  uaer  can  introduce 

data  types  and  operators  between  them.  However,  the  accuracy  and  flexibility  of 

PASCAL-SC  make  it  the  language  of  choice  for  scientific  computation . 

Since  the  implementation  method  described  here  performs  analytic  differentiation 

at  compete,  time,  the  machine  code  produced  is  of  the  same  order  of  efficiency  as  if 

expressions  for  the  derivatives  were  given  in  the  source  code  or,  as  will  be  seen, 

as  if  numerical  differentiation  by  difference  quotients  is  used. 

2.  Mathematical  preliminaries.  The  computation  of  the  value  f t>f (x^,. . . ,xq) 

of  a  function  of  n  real  variables  x, , . . . ,x  is  done  on  a  computer  by  a  sequence  of 

1  n 

arithmetic  operations  and  the  evaluation  of  standard  (or  library)  functions,  for  which 

subroutines  are  available.  From  a  mathematical  point  of  view,  this  means  that  f  is 

the  compoittion  of  a  finite  number  of  functions  f,  ,...,f  ,  that  is, 

X  in 

(2.1)  f  i-  f1«f2»...»fa, 

where  f.  is  in  general  a  function  of  x.,...,x  and  f.  ,,..., f  .  If,  at  the  current 
K  x  n  K+X  fll 

value  of  x:»(x, ,...,x  ),  each  f.  is  differentiable  with  respect  to  its  arguments, 
x  n  k 

then  f  is  differentiable  at  x,  and  f ' (x)  is  obtained  by  the  chain  Ante  of  differential 
calculus , 

(2.2)  f’(x)  s-  f'(x(,"1)).....f  ,  (x{1))*f'(x), 

x  m-x 


where  the  '  denotes  Frtfchet  differentiation  and  the  *  matrix  multiplication  (13 1, 
and 

(2.3)  x(k)  s»  f (k)  (x) ,  f°°  i-  fk+1'...°fm,  *  i-  1 . m-1. 

If  f '  (x)  exists,  then  it  is  represented  by  the  gruuLLent  vectOA 


(2.4)  Vf(x)  -  ...  ,2|i£^T, 

3*1  3*2  3xn 

for  details  about  differentiation  in  vector  spaces,  see  HI],  for  example . 

As  explained  in  (13] ,  once  software  is  in  place  for  the  calculation  of  f^  and 


< 


f£,  then  the  (gradient)  vec-toA  Vf(x)  can  be  obtained  along  with  the  vaJLut  f(x)  at 

x  for  which  f  is  differentiable.  This  software  is  produced  by  the  application  of 

A.U lu  to  the  expressions  defining  f:  The  rules  for  evaluation  of  expressions 

give  the  code  for  f  (x) ,  while  the  rules  fqr  cU^eAWtiation  give  the  corresponding 

code  for  Vf(x).  From  a  practical  standpoint,  all  this  means  is  that  the  result 

of  each  of  the  operations  or  function  evaluations  involved  in  the  computation  of 

f  will  be  the  intermediate  value  of  f,  and  the  corresponding  intermediate  values 

of  the  partial  derivatives  of  f  with  respect  to  x.,...,x  .  Since  the  compiler 

X  n 

is  inforsied  of  the  rules  for  differentiation  of  arithmetic  operations  and  library 
functions,  the  machine  code  is  built  in  the  same  sequential  fashion  as  just  for 
the  evaluation  of  f(x)  only.  The  implementation  of  differentiation  discussed  here 
is  for  a  serial  processor;  as  pointed  out  in  (13) ,  analytic  differentiation  can  be 
done  even  more  efficiently  in  a  parallel  environment. 

To  see  what  will  actually  be  computed  in  mathematical  notation,  suppose  that 
f  is  a  function  of  only  two  independent  variables y,z.  (The  concepts  of  independent 
and  dependent  variables  will  be  defined  precisely  later.)  Then,  if 

(2.5)  f  s«  g(y,z)  , 

where  g (y,z)  denotes  an  expression  containing  only  the  variables  y,z,  then  the  com¬ 
ponents  of  Vf (y ,z)  are  of  course  9f/3y  and  3f/3z  evaluated  at  the  current  values  of 
the  independent  variables.  It  can  happen  that  the  independent  variables  do  not  occur 
explicitly  in  the  expression  defining  f,  for  example,  suppose  that 


(2.6)  f  g(u,v)  , 

where  Vu(y,z)  and  Vv(y,z)  are  known.  Then,  as  usual,  one  gets 

(2  7)  M  -  2SL.UL  +  ,  i£  «  ia.Hi  +  2SL2Z 

3y  3u  3y  3v  3y  3z  3u  3z  3v  3z  ' 

If  one  or  more  independent  variables  also  appear,  such  as  in 
(2.8)  f  g(y,u,v) , 
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then  the  first  component  of  Vf  is 


(2.9) 


ii  »  la  +  ia.lH  +  is.32 

3y  3y  3u  3y  3 v  3y 


an  inaccurate  and  unstable  process,  even  for  functions  of  one  variable,  unless  tech¬ 
niques  of  the  differential  calculus  are  used.  The  user  of  difference  quotients  is 
faced  with  the  traditional  dilessMi  facing  solvers  of  ill-posed  problems:  If  h  is 
too  large,  then  the  difference  quotient  is  a  poor  approximation  to  the  derivative 
(truncation  error),  while  if  h  is  too  small,  then  significant  digits  are  lost  in  the 
maaerator  since  f(x)  and  f(x  he(k))  agree  in  many  of  the  finite  number  of  places 
to  which  they  are  calculated,  because  differentiability  of  f  implies  its  continuity. 
n>us,  as  h  becomes  small,  the  user  of  difference  quotients  is  led  down  the  garden 
path,  losing  one  significant  digit  after  another,  until  one  arrives  in  a  Cloudcukoo- 
land  in  which  the  difference  quotient  is  0  for  all  h  sufficiently  small,  regardless 
of  the  value  of  the  derivative  sought. 

This  situation  creates  a  serious  problem,  which  is  to  find  the  optimal  value 
of  h  in  the  sense  that  as  much  accuracy  as  possible  is  attained  in  the  approximation 


! 


I 


of  the  derivative,  and  the  loae  of  significant  digits  ia  minimized.  An  analysis  by 
Dennis  and  Schnabel  (3]  indicates  the  following  resolution t  If  f(x)  can  be  computed 
with  relative  error  t,  then  h  •  /T  is  a  suitable  choice  in  the  sense  that  the  differ¬ 
ence  quotient  approximates  the  derivative  to  order  t,  while  retaining  about  half  the 
number  of  significant  digits  in  fix),  assuming  the  function  and  the  derivative  being 
computed  are  roughly  the  same  order  of  magnitude,  to  see  how  this  works  out,  it  will 
be  applied  to  the  simple  example 

(3.2)  f  (x)  i»  sinx  , 

a  function  of  a  single  read  variable.  The  method  of  analytic  differentiation  pre¬ 
sented  in  this  paper,  of  course,  gives  the  exact  result 

(3.3)  f'(x)s*lcosx  , 

which  will  be  compared  with  the  difference  quotient 

(3.4)  Ah(f,x)  -hi-:. »-*"■<*> 


for  x  1.00  and  x  j- 1.57.  The  particular  PASCAL-SC  implementation  used  (18)  computes 
with  twelve-digit  decimal  arithmetic;  one  gets 
SIH(I.OO)  -  8.41470984808E-01, 

(3.5) 

SIN{1 .57)  -  9.99999682932E-01. 


Since  these  values  are  guaranteed  to  be  accurate  to  the  twelve  digits  shown  (18] , 

(19),  it  is  reasonable  to  take  t  »  10  12 ,  from  which  the  rule  cited  above  indicates 
-6 

that  h  *  10  is  optimal.  Actual  results  are  shown  in  Tables  3.1  and  3.2. 

The  computed  results  bear  out  the  theory  in  [ 3  ]  to  the  extent  that  what  little 
accuracy  is  obtained  is  maximized  for  the  predicted  value  h  ■  10  The  value  of 
the  exact  derivative,  however,  is  not  only  more  accurate  but  also  computed  fiaiteA 
in  this  case,  even  if  one  also  wants  the  value  of  f(x),  as  is  usually  true.  The  loss 
of  significant  digits  in  fih(f,x)  is  large  in  the  second  case,  where  f(x>  and  f '  (x) 
differ  by  several  orders  of  magnitude,  if  such  a  sacrifice  of  significance  occurs 


-  5  - 


often  enough  in  the  program,  then  the  fine!'  results  will  be  rendered  meaningless 


Higher-order  schemes  for  numerical  differentiation  tic]  suffer  from  the  same  prob¬ 
lems  involving  cancellation  of  significant  digits,  and  are  additionally  slower  to 
compute  than  the  simple  difference  quotient  0.1).  Use  of  analytic  derivatives, 
obtained  automatically,  avoids  this  entire  area  of  problems,  how  that  advanced 
languages  which  permit  easy  implementation  of  differentiation  are  available,  the 
future  of  the  use  of  difference  quotients  to  approximate  derivatives  in  scientific 
computation  appears  to  be  extremely  limited.  It  should  be  kept  in  mind  that  the 
above  remarks  apply  to  differentiation  of  functions  defined  by  computer  programs. 

The  differentiation  of  functions  defined  by  data,  is  a  classical  ill-posed  problem 
[17] ,  for  which  the  mathematical  theory  and  computational  techniques  are  still  being 
developed. 

The  results  in  Tables  3.1  and  3.2  also  show  the  convenience  of  the  accurate 
PASCAL-SC  floating-point  arithmetic  [18] :  the  trailing  zeroB  clearly  indicate  the 
number  of  significant  digits  lost. 

4.  Representation  of  TYPE  GRADIENT.  Returning  to  analytic  differentiation, 
the  discussion  above  shows  that  it  is  natural  to  associate  with  each  differentiable 
function  f  of  n  variables  its  value  f(x)  at  a  point  x,  and  its  n-dimensional  gradient 
vector.  Vf(x)  at  x.  This  association  between  a  real  number  and  a  real  vector  makes 
the  use  of  the  RECORD  data  type  of  PASCAL  appropriate  [ 6 ) .  (PASCAL-SC  is  an  ex¬ 
tension  of  PASCAL,  so  that  concepts  relating  to  the  latter  carry  over  directly.) 

The  standard  declaration  of  TYPE  GRADIENT  then  runs  as  follows: 

CONST  DIM  =  n; 

TYPE  DIMTYPE  =  1 . .DIM; 

<4.1) 

RVECTOR  «  ARRAY  [DIMTYPE]  OF  REAL; 

GRADIENT  *  RECORD  F:  REAL;  DF:  RVECTOR  END; 

The  first  three  lines  of  (4.1)  are  simply  the  declaration  of  TYPE  RVECTOR  in 
PASCAL-SC:  the  value  n  in  lower-case  is  supplied  by  the  user.  Since  n-dimensional 
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4 


f 


i 


real  vectors  are  ubiquitous  in  scientific  computation!  RVECTOR  is  considered  to 
be  a  standard  numerical  data  type  in  PASCAL-SC.  Thus,  by  following  the  form  (4,1) 
exactly,  the  facilities  of  PASCAL-SC  for  computation  with  real  n-dimensional  vectors 
will  be  at  one's  disposal  [18], 

If  F  is  of  type  GRADIENT,  then  F.F  is  called  its  neat  pant,  and  F.DF  its 
vecton  pa/it.  If  F  represents  a  real  function  f  of  x  =  ),  then 


(4.2) 


8f(x) 

3x, 


3f(x) 

3x 


in  the  sense  that  the  corresponding  values  will  be  assigned  to  F.F  and  the  compo¬ 
nents  of  F.DF.  A  more  precise  description  will  be  given  in  the  following  section. 

5.  Representations  of  variables  and  constants.  In  order  to  be  able  to  use 
automatic  differentiation  in  a  sensible  way,  it  is  essential  to  distinguish  between 
independent  and  dependent  variables  on  the  one  hand,  and  constants  on  the  other.  For 
the  purpose  of  differentiation,  all  variables  will  be  of  type  GRADIENT,  while  con¬ 
stants,  as  will  be  seen,  may  be  of  type  INTEGER,  REAL,  or  GRADIENT.  The  necessary 
distinctions  will  be  made  in  more. detail  below. 

5.1.  Independent  variables.  A  variable  V  of  type  GRADIENT  is  said  to 

(K) 

be  the  Kth  tnde.pzn.dwt  vaniable  if  the  Kth  unit  vector  e  is  assigned  to  its 

(K) 

vector  part,  that  is,  V.DF:-e  ,  or,  more  precisely,  if 
(5.1)  V.DF(II:”0  for  I  ?  Kj  V.DF[KJ:»1. 

The  user  of  the  differentiation  software  described  here  is  free  to  name  and  order 
independent  variables  in  an  arbitrary  fashion,  subject  to  the  ordinary  limitations 
of  PASCAL  (6  1.  For  example,  if  x,Y,Z  denote  respectively  the  first,  second,  and 
third  independent  variables,  then 

X. DF(1):=1,  X,  DF  [2]  :=0  ,  X.DF[3]:=0, 

Y. DF[l)!-0,  Y.DF[2]:=1,  Y.DF(3J:=0, 

Z. DF[1J:» 0,  Z.DF(21i=0,  Z.DF[3J:«1, 
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(5.2) 


define  their  vector  parts  X.DF,  Y.DF,  and  2. DP,  respectively.  If  independent  var¬ 
iables  appear  in  assignment  statements,  it  should  be  only  on  the  right  side.  An 
assignment  to  the  vector  part  of  the  Kth  independent  variable  of  anything  but  the 
Kth  wit  vector  makes  that  variable  dependent  (see  below) ,  and  introduces  a  new, 
but  perhaps  unspecified  Kth  independent  variable,  with  respect  to  which  subsequent 
partial  derivatives  will  be  computed.  Since  the  rules  for  partial  differentiation 
are  applied  rapidly  and  accurately,  but  mindlessly,  during  the  course  of  the  compil¬ 
ation,  the  user  should  avoid  introducing  errors  of  this  type  wless  a  change  of 
variable  is  actually  intended.  In  the  latter  case,  assignments  to  former  independent 
variables  of  their  values  in  terms  of  new  ones  performs  the  change  of  variables 
automatically,  given  the  corresponding  expressions  and  subroutines.  Of  course,  it 
will  often  happen  in  the  course  of  the  computation  that  assignments  will  be  made  to 
the  value  part  V.F  of  an  independent  variable  V. 

Another  possible  source  of  difficulty  in  setting  up  independent  variables  would 
be  the  assignment  of  the  same  wit  vector  to  the  vector  parts  of  two  supposedly  in¬ 
dependent  variables,  say  V  and  W,  that  is,  V.DF  -  W.DF.  This  mews  that  V  and  W 
would  be  considered  to  be  the  same  variable  for  the  purpose  of  differentiation,  but 
with  possibly  different  values  V.F  and  H.F.  This  might  be  used  as  a  sneaky  way  of 
assigning  several  values  to  the  Kth  independent  variable,  but  is  not  recommended, 
since  the  same  purpose  can  be  served  by  a  straightforward  method  which  avoids  the 
possibility  of  confusion  and  perhaps  well -concealed  errors. 

5.2.  Dependent  variables.  Variables  appearing  on  the  left  side  of 
assignment  statements  are,  of  course,  dzpe.nde.rtf  on  the  variables  appearing  in  the 
expressions  on  the  right  side,  and  thus  ultimately  on  the  set  of  independent  var¬ 
iables  chosen  by  the  user.  A  dependent  variable  F  depends  on  the  Kth  independent 
variable  or  not  according  as  F.DF [K)  ^  0  or  F.DFfK)  »  0. 

5.3.  constants.  A  "variable"  of  type  GRADIENT  is  said  to  be  a  c on&tant 
if  its  vector  part  is  the  n-dimensional  zero  vector  0  «  (0,...,0),  that  is,  C  is  a 
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constant  if 


(5.3)  C.DF [I]  -  0,  I  -  l,...,n. 

Constants  of  type  GRADIENT  can  be  produced  in  the  computation  by  direct  assignment, 
or  by  expressions  such  as  F:-0*X,  F:»x/x,  F:«X-X,  and  so  on.  It  is  also  permissible 
to  have  variables  or  literal  constants  of  type  INTEGER  or  REAL  in  GRADIENT  expressions; 
such  quantities  will  be  treated  as  constants  for  the  purpose  of  differentiation. 

Thus,  an  assignment  statement  such  as 

(5.4)  F  i-  <K  +  X)  **(4.4  +  V)  -  R  -  3; 

where  K  is  of  type  INTEGER,  R  is  REAL,  and  F,X,Y  are  GRADIENT,  is  allowed.  This 
makes  it  easy  to  handle  expressions  involving  real  or  integer-valued  parameters. 

If  X  and  Y  are  independent  variables,  then  only  3F/3X  and  3F/3Y  will  be  altered  in 
F.DF  by  the  assignment  (5.4),  in  addition  to  the  value  part  F.F  of  F. 

6.  Definition  of  operators  for  type  GRADIENT.  In  order  to  compute  with 
variables  of  type  GRADIENT,  the  arithmetic  operators  and  the  standard  and  other 
functions  needed  must  be  defined  in  order  to  conform  to  the  rules  for  evaluation 
and  differentiation  of  real  functions.  Arithmetic  operations  +,-,*,/,*•  will  be 
considered  in  this  section,  and  the  standard  and  other  functions  in  the  next.  In 
what  follows,  x,  R,  ra,  RB,  G,  GA,GB  will  denote  generic  variables  of  the  follow¬ 
ing  types: 

(6.1)  VAR  X:  INTEGER;  R.RA.RB:  REAL;  G,GA,0:  GRADIENT; 

the  operator  extension  in  PASCAL-SC  [  1  ]  will  be  employed  to  obtain  the  desired 
results.  The  source  code  for  the  arithmetic  operators  is  given  in  Appendix  A. 

6.1.  Addition  (+) ■  A  total  of  six  definitions  of  the  addition  opera¬ 
tors  are  required:  One  each  for  the  various  combinations 

(6.2)  +G,  X+G,  G+X,  R+G,  0+R,  GA+GB. 

For  example,  the  unary  addition  (identity)  operator  denoted  by  +G  is  declared  by 
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♦  / 


(6.3) 


OPERATOR  *  (G:  GRADIENT)  RES:  GRADIENT; 


BEGIN  RES:=G  END; 

according  to  the  rules  of  PASCAL-SC  (1).  Similarly,  the  assignments  F i "K+G ,  F:«G+K, 
F : “R+G ,  F : *G+R  alter  only  the  value  part  of  Fj  one  has  respectively  F.F««K+G.F, 
F.F:»G.F+K,  F.Fi“fttG.F,  F.F:“G+R,  while  F.DFi^G.DF  in  each  case.  Finally,  corres¬ 
ponding  to  GA+GB,  one  has  the  operator  declaration 


OPERATOR  +  (GA:  GRADIENT ;GB:  GRADIENT)  RES:  GRADIENT; 
VAR  U:  GRADIENT;!:  DIMTYPE; 

BEGIN  U.F:=GA.F+GB.F;FOR  I:-l  TO  DIM  DO 

(6.4) 

U.OF[I]:=GA.DF[I]+GB.DF[I]; 

RES:*U 

END; 


since  the  derivative  of  the  sum  is  the  sum  of  the  derivatives.  In  mathematical 
notation,  if  u,v  denote  generic  functions,  c  a  constant,  and  x  an  independent  var¬ 
iable,  the  operator  declarations  simply  implement  the  rules  for  differentiation 


(6.5) 


£(+u)  .  ±{c  ♦  u)  -  £<u  +  c) 


3u 
3x  ' 


^tlU+V) 


3u  3v 
3x  3x 


for  the  vector  parts  of  dependent  gradient  variables. 

6.2.  Subtraction  (-) .  Once  again,  operators  are  required  for  the  com¬ 
binations 


(6.6)  -G,  K-G,  G-K,  R-G,  G-R,  GA-GB 

of  GRADIENT,  INTEGER,  and  READ  types.  For  example,  unary  subtraction  -G  (sign 
changing) ,  is  accomplished  in  GRADIENT  statements  by  the  operator  declared  by 


(6.7) 


OPERATOR  -  (G:  GRADIENT)  RES:  GRADIENT; 

VAR  U:  GRADIENT;!:  DIMTYPE; 

BEGIN  U.F:*-G.F;F0R  I:=l  TO  DIM  00  U.DF[IJ:*-G.DF[I];RES:*U 
END; 
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Here,  the  rules  corresponding  to  (6.5)  are 


(6.8) 


3u 

3x' 


3u 

3x' 


3u  3v 
3x  “  3x  ' 


which,  in  symbolic  fora,  give  the  assignments 


(K-G) .Fi«K-G.F>  (K-G)  .DFi— G.DF;  (R-G) .F:-R-G.F>  (R-G) .DF:»- G.DF, 

(6.9)  (G-K).Fs-G.F-K*  (G-K) .DFi-G.DFj  (G-R) .Fi-G.F-R,  (G-R) .DF:«G.DF, 

(GA-GB) ,F:-GA. F-GB.F;  (GA-GB) . DF : “GA.DF-GB.DF  j 


where  the  vector  operations  are  understood  to  be  performed  componentwise,  as  in 
usual  vector  algebra. 

6.3.  Multiplication  (*) .  Here,  operators  must  be  defined  for  the  com¬ 
binations 


(6.10)  K*G,  G*K,  R*G,  G*R,  GA*GB. 


The  rules  for  differentiation  of  products  are, 
(6.11)  £<c*u)  -  gtu-c)  -  eg.  g(u.v) 


of  course. 


3v  ^  3u 
U*  3x  3x  v' 


which  lead  to  the  symbolic  assignments 

(K*G)  .F:«K*G.F,  (K*G) .DFi-K*G.DF»  (R*G) .F:-R*G.F»  (R*G) ,DF:-R*G.DFi 

(6.12)  (G«K) .F;-G.F«K,  (G*K) ,DF:-G.DF*X|  (G*R) .F:“G.F*R|  (G*R)  .DF:-G.DF»R, 

(GA*GB) .F:«GA.F*GB.F;  (GA*GB) .DF:-GA.F»GB.DF+GA.DF*GB.F» 


the  detailed  code  is  given  in  Appendix  A. 

6.4.  Division  (/) ■  Here,  as  for  multiplication,  the  combinations  for 
which  operators  are  needed  are 


(6.13)  K/G,  G/K,  R/G,  G/R,  GA/GB. 

The  rules  for  differentiation  of  quotients  are  implemented  in  the  forms 

(6.14)  g(c/u>  -  -<c/u)-g/u,  g(u/c)  -  g/c,  g(u/v)  -  g/v  -  g-(u/v)/v. 
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I 


Letting  C  stand  for  K  or  R,  these  correspond  to  the  assignments 


(C/G).F:-C/G.F>  (C/G)  .DF G.DF*  (C/G)  .F/G.F; 

(6. 15)  (G/C) .F:-G.F/C|  (G/C) .DF:-G.DF/K) 

(GA/GB) .Fs-GA.F/GB.F)  (GA/GB) ,DF:“GA.DF/GB.F-GB.DF* (GA/GB) . F/GB.Fj 


in  symbolic  form.  Explicitly,  source  code  for  the  operator  /  between  gradient  op¬ 
erands  GA,GB  is: 


(6.16) 


OPERATOR  /  (GA:  GRADIENT ;GB:  GRADIENT)  RES:  GRADIENT; 

VAR  U:  GRADIENT ; I :  DIMTYPE; 

BEGIN  U.F:«GA.F/GB.F;FOR  I:*1  TO  DIM  DO 

BEGIN  IF  GA.DF[I]  =■  0  THEN  U.DF[I]:«0  ELSE  U.DF[I] :=GA.DF[I]/GB.F 
IF  GB.DFCI]  <>0  THEN  U.DF[I] :«U.DF[I]-GB.DF[I]*U.F/GB.F 

END; 

RES:*U 

END: 


as  given  in  Appendix  A.  Attempted  division  by  zero  produces  an  error  interrupt. 

6.5.  Power  (**) .  The  power  operator  defined  by 

(6.17)  0**v  -  uV 

is  not  standard  in  PASCAL  or  PASCAL-SC.  Therefore,  this  operator  is  also  introduced 
for  R**K,  a  REAL  raised  to  an  INTEGER  power,  and  RA* *RB  for  RA.RB  of  type  REAL.  An 
algorithm  which  is  recommended  for  this  purpose  [  2 )  is  based  on 


(6.18) 


U*»V 
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V*LOG10(U) 


since  decimal  arithmetic  is  being  used.  However,  experiments  with  integer  exponents 
showed  that  more  accuracy  was  obtained  in  this  case  by  use  of  the  method  from  [13] 
of  repeated  squaring  of  the  base,  as  employed  in  the  original  differentiation  soft¬ 
ware  written  by  Reiter  [  5  }.  Specifically,  assume  that  K  is  a  non-negative 
integer,  and 
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(6.19) 


K  -  eQ  +  6^2  ♦  e2.22  ♦  ...  +  V2B, 
where  £_  »  1  and  e .  ■  0  or  1,  i  -  0,l,...,m-l.  Then, 

K  “  e  -21  “  E 

(6.20)  R  -  n  Ri  «  n  R.  i, 

i-0  i-0  1 

where 

2 

(6.21)  Rq  «  R,  Ri  “  i*  i  ”  l,...,m. 

Thus,  each  factor  Ri  is  formed  by  squaring  the  previous  R^,  and  is  multiplied  into 

the  product  (6.20)  or  not  according  «i  e.  *  1  or  0.  Hence  the  name  "repeated  squaring' 

for  this  algorithm.  For  negative  K,  the  reciprocal  of  the  result  for  positive  K  is 

0  K 

taken.  As  pointed  out  in  (13) ,  one  takes  R  »  1  for  R  ?  0,  0  -  0  for  K  >  0,  while 

K 

0  for  K  s  o  leads  to  an  error  interrupt  due  to  an  attempted  division  by  0.  The  case 
Rl  ■  R  is  also  handled  separately  in  the  source  code,  which  is 
OPERATOR  **  (R:  REAL;K:  INTEGER)  RES:  REAL; 

VAR  L:  INTEGER ;U:  REAL; 

BEGIN  IF  K  (»  0  THEN  U:*1/R;IF  K  =  0  THEN  U:«l 
ELSE  IF  K  *  1  THEN  U:*R 

ELSE  BEGIN  L:=ABS(K) ;U:*1 ;REPEAT  IF  L  MOD  2  *  1 

(6.22)  THEN  0:=R*U;L:*L  DIV  2;  IF  L  <>  0 

THEN  R:»R*R  UNTIL  L  *  0; 

IF  K  <0  THEN  U:»l/U 

END; 

RES :  *1) 

END; 

The  computation  of  RA**RB  for  a  REAL  exponent  RB  is  based  on  the  same  algorithm;  one 
computes  RA»*K  for  the  INTEGER  exponent  K:-TRUNC(RB) .  Then,  -1  <  RB  -  K  <  1,  and 
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<6.23) 


RA*»RB:-(RA«*K) «EXP ( (RB-K) *LN (RA) ) ; 


$ 


t 


in  which  natural  logarithms  are  used.  Negative  bases  RA  with  nonintegral  exponents 
will  lead  to  an  error  interrupt,  since  the  logarithm  function  will  not  accept  neg- 

DO 

ative  arguments i  an  attempt  to  compute  0  for  RB  s  0  will  similarly  be  wrecked  by 
a  division  by  zero  error  interrupt. 

Once  operators  for  integral  and  real  powers  of  real  arguments  are  in  place, 
the  corresponding  operators  for  gradient  variables  can  be  derived  from  the  rules 
for  differentiation , 


(6.24) 


3  / 

3?(u  > 


c-1 


1»  _L(eu) 

3x'  3xl  ' 


u  ,  ,  .  3u 
c  *ln  (c)  , 


and 


(6.25) 


3  ,  v, 
3l<U  > 


v-1  3u  ,  v  ,  ,  ,  3v 
v>u  •—  +  u  »ln(v)*— . 
3x  3x 


The  power  operator  ••  is  thus  available  for  the  combinations 

(6.26)  R**K,  RA**RB,  K*»G,  G*»X,  R**G,  G**R,  GA**GB. 

For  the  forms  involving  GRADIENT  expressions,  the  symbolic  assignments  are t 

(G»*C) ,F:-G.F**C»  <G**C) .DF:-C»( (G»*C) .F/G.F) *G.DF; 

(6.27)  (C**G) .F:-C**G.F>  (C**G) .DF«-(C»*G.F) *LN(C) *G.DF, 

(GA**GB) .F:»GA.F*«GB.Fj  (GA**GB) .DFs“(GA.F**GB) . DF+ (GA*  *GB . F) . DF; 

once  again,  all  operations  involving  the  vector  parts  .DF  of  the  gradient  variables 
are  to  be  interpreted  componentwise. 

6.6.  Priorities  of  Operators.  The  introduction  of  type  GRADIENT  thus  re¬ 
quires  the  definition  of  29  operators;  six  each  for  addition  and  subtraction,  five 
each  for  multiplication  and  division,  and  the  seven  power  operators,  including  the 
two  for  raising  a  READ  base  to  an  INTEGER  or  REAL  power.  The  latter  can  also  be 
used  in  ordinary  numerical  computation.  The  priorities  of  these  gradient  operators 
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and  the  two  real  operators  defined  above  are  as  follows: 


1st  priority:  unary  addition  and  subtraction  ±; 

2nd  priority:  Multiplication,  division,  and  power  *,/,**» 

3rd  priority:  Binary  addition  and  subtraction  +,-. 

As  usual,  parentheses  are  introduced  to  achieve  the  desired  order  of  operations. 

Users  familiar  with  languages  in  which  **  has  a  higher  priority  than  *,/  should  be 
particularly  careful  to  make  their  desires  explicit.  For  example,  2*3**4  Mans 
64  -  1,296,  not  2-34  -  162. 

7.  Standard  gradient  functions.  The  present  implementation  for  type  GRADIENT 
includes  the  following  standard  function:  Absolute  value,  and  the  six  standard  func¬ 
tions  for  type  REAL  available  in  the  PASCAL-SC  compiler  [18)  (square  root,  exponential 
(base  e) ,  natural  logarithm,  arctangent,  sine,  and  cosine).  All  standard  gradient 
factions  have  names  which  begin  with  G:  GABS,  GSQRT,  GEXP,  GIN,  GARCTAN ,  GSIN, 
and  GOOS.  Thus,  to  obtain  the  value  and  the  gradient  vector  for  the  function 

(7.1)  f(x.y)  -  (xy  +sinx  +  4)  (3y2  +  6), 

til),  [12 1  »  (13)  ,  the  corresponding  GRADIENT  assignment  statement  is 

(7.2)  F:-(X*Y  ♦  GSIN(X)  +■  4>*(3«(Y»»2)  ♦  6); 

keeping  in  mind  the  equal  priority  of  *  and  **.  (The  software  described  in  this 
paper  was  tested  first  on  (7.2),  for  historical  reasons.)  The  minor  nuisance  of 
having  to  write  GSIN(X)  instead  of  SIN(X)  is  balanced  by  the  fact  that  the  forssr 
makes  it  clear  that  (7.2)  is  a  GRADIENT  assignment  stateMnt,  and  GRADIENT  arguMnts 
are  expected  in  the  expression  on  the  right,  at  least  for  X  in  the  sine  function 
and  thus  also  in  the  product  X*Y. 

For  the  absolute  value,  one  has 

(7.3)  ^|u|  ■  -  if  u  <  0,  gjj:|u|  -  |j-  if  u  >  0,  ^j|u|  undefined  if  u  -  0 

In  source  code  (see  Appendix  B) ,  this  becomes : 
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FUNCTION  GA8S(G:  GRADIENT):  GRADIENT ; 


VAR  I:  DIHTYPEjN:  REAL ;U :  GRADIENT; 

BEGIN  U.F:*AB$(G.F);M:*G.F/U.F;FOR  I:-  1  TO  DIM  DO  IF  G.DF[I]«0  THEN 

(7.4) 

U.DF[I]:«0  ELSE  U.DF[I]:«M*G.DF[I]; 

GABS:-U 

END; 

the  nondiffarentiablity  of  this  function  at  zero  will  be  indicated,  if  attempted, 
by  a  division  by  sero  error  interrupt. 

The  remaining  standard  functions  follow  the  pattern  (7.4)  exactly,  using  their 
known  derivatives  to  construct  assignments  to  the  mittiptieA.  (or  divisor)  M.  In 
brief  form,  the  required  assignments  are 

U.Fi-SQRT(G.F);  M:-2*U.F»  U.DF!-G.DF/M»  GSQRT:-U> 

U.Fj-EXP(G.F);  U.DF:«U.F*G.DF»  GEXPs-U; 

U.Fi-UI(G.F)»  U.DF.--G.DF/G.F;  GLN:«Uf 

(7.5) 

O.F«-ARCTAN(G.F);  Ms-1-MS.F*G.F»  O.DF.--G.DF/M;  GARCTAN:=*U» 

U.Ft-SIN(G.F)  *  M:-OOS(G.F)f  U.DF.—M*G.DF»  GSIN :“U; 

U.F:*OOS (G.F)  ;  M:— SIN(G.F);  U.DF«-M«G.DF;  GOOS.—U; 

for  the  respective  standard  functions.  An  explicit  assignment  to  M  was  unnecessary 
in  the  case  of  the  exponential  and  logarithmic  functions. 

8.  User-defined  functions.  The  user,  of  course,  is  free  to  add  functions 
(or  procedures)  for  type  GRADIENT,  as  long  as  the  rules  for  differentiation  of  the 
results  are  known  explicitly.  In  the  case  of  functions  of  a  single  variable,  the 
code  (7.4)  can  act  as  a  template  for  the  required  code.  If  the  real  function  f 
is  known  as  FUNC,  and  its  derivative  f'  is  computed  by  DFUNC,  then  the  gradient 
function  GFUNC  requires  the  assignments 

(8.1)  U.Fi-FUNC(G.F) i  M: -DFUNC (G. F) >  U.DF:-M*G.DF>  GFUNCi«U» 


which  correspond  to  the  rule  for  partial  differentiation 

(8.2)  ^f(u)  -  f  (u)*|j. 

For  example ,  suppose  that  the  cosecant  and  its  derivatives  are  needed.  Then,  fol¬ 
lowing  the  pattern  of  (7.4),  the  user  can  introduce  the  function  GCSC  by  the  code: 

FUNCTION  GCSC(G:  GRADIENT):  GRADIENT; 

VAR  I:  DIMTYPE ;M:  REAL ;U :  GRADIENT; 

BEGIN  U.F:-1/S1N(G.F);M:--C0S(G.F)*U.F*U.F;F0R  I:-l  TO  DIM  DO 

(8.3) 

IF  G.DF[I]  -  9  THEN  U.DF[I]:-0  ELSE  U.DF[I]:*M*G.DF[I]; 

GCSC:-U 

END; 

however,  since  the  function  in  question  can  be  expressed  easily  in  terms  of  available 
gradient  functions  and  operators,  a  more  compact  code  is 

FUNCTION  GCSC(G:  GRADIENT):  GRADIENT; 

(8.4) 

BEGIN  GCSC:«1/GSIN(G)  END; 


The  coding  of  gradient  functions  and  operators  can  also  be  simplified  by 
using  the  operators  and  functions  for  vector  and  matrix  algebra  available  in  the 
PASCAL-SC  library  [18] .  This  approach  is  particularly  useful  if  it  is  desired  to 
precompile  the  gradient  operators  and  functions ,  and  store  the  relocatable  code  gen¬ 
erated  in  an  external  library.  This  cannot  be  done  for  the  source  code  given  in 
the  appendices,  since  most  of  these  routines  contain  the  global  variable  DIM,  which 
can  be  removed  if  source  code  for  vector  operations  is  brought  into  the  program. 

In  vector  form,  (8.3)  becomes 

FUNCTION  GCSC(G:  GRADIENT):  GRADIENT; 

VAR  M:  REAL;U:  GRADIENT; 

(8.5)  BEGIN  U.F:-1/SIN(G.F);M--C0S(G.F)*U.F*U.F;U.DF:-M*G.DF; 

GCSC:*U 

ENO; 
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with  no  global  variables.  In  all  cases,  the  assignment  U.DF:“M*G.DF  should  be 
arranged  so  that  the  multiplication  of  the  RVECTOR  G.DF  by  the  REAL  M  is  frost  the 
left. 

9.  An  example!  newton's  method  for  the  solution  of  nonlinear  systems  of 
equations.  Type  GRADIENT  has  many  applications  in  scientific  computing,  for  exam¬ 
ple,  to  sensitivity  and  error  analysis,  optimisation,  and  the  solution  of  systems 
of  nonlinear  equations  (13] .  In  order  to  apply  Newton's  method  to  the  numerical 
solution  of  the  system  of  equations, 

£l(xl’x2""'xn>  “ 
fj (w^ **2 • *  *  * 'Xn^  * 

(9.1) 


f  (x  ,x  ,  • » « ,x  )  *  0, 
n  l  2  n 


the  n»n  Jacobian  matJvix  J  “  (3fi/3x^)  is  needed  (11] .  Without  automatic  differen¬ 
tiation  ,  the  labor  involved  in  producing  J  from  (9.1),  not  to  speak  of  the  possibility 
of  introduction  of  errors,  is  prohibitive  even  for  moderate  n  in  the  general  case. 
However,  if  Fl,...,FN  and  Xl,...,XN  are  of  type  GRADIENT,  then  the  Ith  row  of  the 
Jacobian  matrix  J  is  simply  FI.DF,  the  vector  part  of  FI,  assuming  Xl,...,XN  are 
independent  variables.  Specifically,  suppose  that  the  PASCAL-SC  type  RMATRIX  is 
declared  by 

(9.2)  TYPE  RMATRIX  =  ARRAY  [DIMTYPE]  OF  RVECTOR; 

where  DIMTYPE  and  RVECTOR  are  introduced  as  in  (4.1),  and  the  new  type  JACOBIAN  by 

(9.3)  TYPE  JACOBIAN  «  ARRAY  [DIMTYPE]  OF  GRADIENT; 

further,  suppose  X,F  are  of  type  Jacobian,  where  X(I] .DF  is  the  Ith  unit  vector. 

Then,  X(I]  will  be  the  Ith  independent  variable,  and  for 

(9.4)  FlI]:-fI(X(l],X(2],...X(n]),  I«l,....n, 

the  problem  of  solving  the  system  (9.1)  amounts  to  finding  the  solution  vector  with 
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components  x(I  J  .F, . . .  ,X  tn]  .F  such  that  F[I].F  -  0,  I  -  l,...,n. 

T  T 

Returning  to  the  system  (9.1),  let  x  “  ,  f(x)  ■  (f^(x),...,  £^(x) )  . 

Then,  one  step  of  Newton's  method,  starting  from  the  approximation  x*°*  to  the  solu¬ 
tion  x  of  (9.1)  ,  requires  the  solution  of  the  linear  system  of  equations 

(9.5)  J*«  -  -f(x<0)) 

T  (1) 

for  the  c oKXiction  vector  5  -  («1,..., &n)  ,  from  which  the  next  approximation  x 

to  x  is  obtained  as 


see,  for  example ,  (11].  Suppose  that  one  has  a  function 

(9.7)  FUNCTION  SOLVLN  (N:  DIMTYPEjA:  RMATRIX;B:  RVECTOR):  RVECTOR; 


which  will  give  the  solution  of  a  system  of  N  linear  equations  in  N  unknowns  with 
coefficient  matrix  A  and  right-hand  side  B.  (A  very  accurate  procedure  for  the 
solution  of  linear  systems  called  LGLP  is  available  in  PASCAL-SC  [18 1.)  The  code 
for  one  step  of  Newton's  method  with  X,F  of  type  JACOBIAN  and B ,D  of  type  RVECTOR, 

J  (the  Jacobian  matrix) ,  of  type  RMATRIX,  is 

FOR  I :*1  TO  DIM  DO  BEGIN  j[l]  :« F[I].DF;B[I]:«-F[I].F  END; 

(9.8)  D:=S0LVLN(DIM,J,B) ; 

FOR  I:*I  TO  DIM  DO  X[I].F:=XtI].F+D[I]; 

for  the  general  case. 

It  is  not  necessary  or  always  convenient  to  introduce  the  type  JACOBIAN. 

For  example,  suppose  that  X,y,Z  are  the  independent  variables  of  type  GRADIENT  (see 
(S.2) ,  and  the  GRADIENT  dependent  variables  F,G,H  are  defined  by 

F:“16* (X**4)+16* (Y**4) +Z**4-16 ; 

(9.9)  G:-X»*2+Y**2+Z**2-3; 

H:»X**3-Y» 


thus,  the  conditions  F.F  -  G.F  =  H.F  "  0  correspond  to  the  system  of  equations 
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* 


(9.10) 


i 

i 


16  x^  +  16y4  +  t*  -  16  -  0, 
x2  +  y2  +  *2  -  3-0, 
x3  -  y  -  0, 


which  waa  chosen,  again  for  hiitorical  reasons  [4],  to  test  the  gradient  Newton 
method  program.  For  the  assignments  (9.9) ,  the  code  (9.8)  is  expanded  for  B,D  of 
type  KVECTOR,  J  (the  Jacobian  matrix)  of  type  RMATRIX ,  to 

J[1]:»F.DF;J[23:«G.DF;J[3]:-H.DF; 

B[1]:-F.F;B[2]:-6.F;B[3]:--H.F; 

(9.11) 

D:-S0LVLN{DIM,J,B); 

X.F:»X.F+D[l]jY.F:«Y.F+D[2];Z.F:«Z.F+D[3]; 


the  result  of  eight  iterations,  starting  from  X.F  -  Y.F  -  Z.F  -  1,  is 

X.F  -  8 . 77965760274E-01 ,  F.F  -  0 . OOOOOOOOOOOE+OO , 

(9.12)  Y.F  -  6.76756970518E-01 ,  G.F  -  0 . OOOOOOOOOOOE+OO , 

Z.F  -  1 . 3308SS41X62E+00 ,  8.F  -  0. OOOOOOOOOOOE+OO. 

a  iiiimsty  of  the  complete  calculation  is  given  in  appendix  C.  the  results  agree 
exactly  to  the  given  number  of  places  to  the  results  for  (9.10)  obtained  on  a 
UN I vac  1100  in  double  precision,  using  the  program  NEKTON  [ 8 ] ,  which  also  employs 
automatic  differentiation. 

10.  implenentatlon  details.  The  software  described  in  this  report  was  im¬ 
plemented  and  tested  on  a  Zilog  MCZ  1/05  microcomputer.  This  machine  has  a  Z80 
processor  and  64  kilobytes  of  main  storage,  augmented  by  two  8"  disk  drives  for 
single-sided,  single-density  hard-sectored  disks.  With  one  drive  dedicated  to  the 
system  disk,  this  gives  the  user  about  308  kilobytes  of  mass  storage.  The  PASCAL-SC 
compiler  used  was  obtained  from  the  Institute  for  Applied  Mathematics,  university 
of  Karlsruhe,  Germany,  and  rims  under  the  Zilog  RIO  2.06  operating  system.  The  modest 
sise  of  the  machine  limits  one  to  about  20«20  matrices.  The  program  to  solve  the 
system  (9.10)  required  less  than  5,632  bytes,  showing  the  compactness  of  the  software. 
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XI.  Conclusions.  By  an  actual  example ,  it  has  been  shown  that  a  modem 
computer  language  which  permits  user-defined  operators  and  data  types  allows  the 
implementation  of  automatic  analytic  differentiation  to  compute  the  gradient  vector 
as  well  as  the  value  of  real  functions  of  n  real  variables .  This  allows  the  accuracy 
and  theoretical  advantages  of  the  use  of  analytically  defined  derivatives  in  actual 
computation  [12 ]<  (13 1.  and  avoids  the  problems  associated  with  the  approximation 
of  derivatives  by  difference  quotients  {  3  ] .  In  addition ,  source  code  using  type 
GRADIENT  is  easier  to  prepare  and  understand  than  if  approximations  to  derivatives 
are  used,  and  the  compiled  code  is  comparable  in  size  and  execution  time  to  that 
produced  in  the  latter  case. 
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APPENDIX  A.  SOURCE  CODE  FOR  GRADIENT  OPERATORS 


(Pile  GRAD  ARITH.S) 


A.l.  Addition  Operators  (Pile  GRAD_ADD.S) 

OPERATOR  +  (G:  GRAOIENT)  RES:  GRADIENT; 
BEGIN  RES:=G  END; 


OPERATOR  +  (K:  INTEGER ;G:  GRADIENT)  RES:  GRADIENT; 
VAR  U:  GRADIENT; 

BEGIN  U.F:=K+G.F;U.DF:=G.DF;RES:*U  END; 


OPERATOR  +  (G:  GRADIENT ;K:  INTEGER)  RES:  GRADIENT; 
VAR  U:  GRADIENT; 

BEGIN  U . F : *G. F+K ;U . DF : «G . DF : RES : =U  END; 


OPERATOR  +  (R:  REAL;G:  GRADIENT)  RES:  GRADIENT; 
VAR  U:  GRADIENT; 

BEGIN  U.F:-R+G.F;U.DF:=G.DF;RES:=U  END; 


OPERATOR  +  (G:  GRADIENTS:  REAL)  RES:  GRADIENT; 
VAR  U:  GRADIENT; 

BEGIN  U.F:«G.F+R;(J.OF:»G.DF;RES:-U  END; 


OPERATOR  +  (GA:  GRADIENT ;GB:  GRADIENT)  RES:  GRADIENT; 
VAR  U:  GRADIENT ; I :  DIMTYPE; 

BEGIN  U.F:=GA.F+GB.F;FOR  I:«1  TO  DIM  DO 
U . DF[ I ] : =GA . DF[ I ]+GB . DF [ I ] ; 

RES:*U 

END; 


A. 2.  Subtraction  Operators  (File  GRADjSUB.S) 

OPERATOR  -  (G:  GRADIENT)  RES:  GRADIENT; 

VAR  U:  GRADIENT; I:  DIMTYPE; 

BEGIN  U.F:--G.F;FOR  I:»l  TO  DIM  DO  U.DF[I]:»-G.DF[I];RES:-U 
END; 
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OPERATOR  -  (K:  INTEGER*:  GRADIENT)  RES:  GRADIENT; 

VAR  U:  GRADIENT ; I :  DIMTYPE; 

BEGIN  U.F:«K-G.F;FOR  I :-l  TO  DIM  DO  U.DF[I]:--G.DF[I];RES:«U  END; 

OPERATOR  -  (G:  GRADIENT*:  INTEGER)  RES:  GRADIENT; 

VAR  U:  GRADIENT; 

BEGIN  U.F:*G.F-K;U.DF:«G.DF;RES:*U  END; 


OPERATOR  -  (R:  REAL;G:  GRADIENT)  RES:  GRADIENT; 

VAR  U:  GRADIENT ;I :  DIMTYPE; 

BEGIN  U.F:*R-G.F:FOR  I:»1  TO  DIM  DO  U.DF[I]:«-G.DF[I];RES:«U  END; 

OPERATOR  -  (G:  GRADIENT*:  REAL)  RES:  GRADIENT; 

VAR  U:  GRADIENT ; 

BEGIN  U.F:-G.F-R;U.DF:-G.OF;RES:*U  END; 

OPERATOR  -  (GA:  GRADIENT ;GB:  GRADIENT)  RES:  GRADIENT; 

VAR  U:  GRADIENT ;I :  DIMTYPE; 

BEGIN  l).F:«GA.F-GB.F;FOR  I:-1  TO  DIM  DO 
U.DF[I]:-GA.DF[I]-GB.OF[I]; 

RES:*U; 

END; 

A. 3.  Multiplication  Operators  (File  GRAD_MUL.S) 

OPERATOR  *  (K:  INTEGER*:  GRADIENT)  RES:  GRADIENT; 

VAR  U:  GRADIENT; I:  DIMTYPE; 

BEGIN  U.F:«K*G.F;FOR  I:-  1  TO  DIM  DO  U.DF[I]:-K*G.DF[I];RES:«U  END 

OPERATOR  *  (G:  GRADIENT*:  INTEGER)  RES:  GRADIENT; 

VAR  U:  GRADIENT ;I:  DIMTYPE; 

BEGIN  O.F:-G.F*K;FOR  I:«l  TO  DIM  DO  U.DF[I]:«G.DF[I]*K;RES:»U  END; 

OPERATOR  *  (R:  REAL*:  GRADIENT)  RES:  GRADIENT; 

VAR  U:  GRADIENT;!:  DIMTYPE; 

BEGIN  U.F:«R*G.F;FOR  I:«l  TO  DIM  DO  U.DF[I]:-R*G.DF[I];RES:«U  END; 
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OPERATOR  *  (G:  GRADIENT-.R:  REAL)  RES:  GRADIENT; 

VAR  U:  GRADIENT ;I:  DIMTYPE; 

BEGIN  U.F:«G.F*R;FOR  I:-l  TO  DIM  DO  U.DF[I]:«G.DF[I]*R;RES:«U  END; 

OPERATOR  *  (GA:  GRADIENT ;GB:  GRADIENT)  RES:  GRADIENT; 

VAR  U:  GRADIENT;!:  DIMTYPE; 

BEGIN  U.F:-GA.F*GB.F;FOR  I:-l  TO  DIM  DO 

U.DF[I]:-GA.F*GB.DF[I]+GA.DF[1]*6B.F; 

RES:*U 

END; 


A. 4.  Division  Operator*  (pile  GRAD_DIV.S) 

OPERATOR  /  (K:  INTEGER ;G:  GRADIENT)  RES:  GRADIENT; 

VAR  U:  GRADIENTS:  DIMTYPE; 

BEGIN  U.F:*K/G.F;FOR  I :-l  TO  DIM  DO 

IF  G.DF[I]  -  0  THEN  U.DF[I]:-0  ELSE  U.DF[I]:«-G.DF[I]*U.F/G.F; 
RES:-U 

END; 

OPERATOR  /  (G:  GRADIENT ;K:  INTEGER)  RES:  GRADIENT; 

VAR  U:  GRADIENT;!:  DIMTYPE; 

BEGIN  U.F:*G.F/K;FOR  I:-1  TO  DIM  DO  U.DF[I]:»G.DF[I]/K;RES.-U  END; 

OPERATOR  /  (R:  REAL ;G :  GRADIENT)  RES:  GRADIENT; 

VAR  U:  GRADIENT;!:  DIMTYPE; 

BEGIN  U.F:*R/G.F;FOR  I:-1  TO  DIM  DO 

IF  G.DFCI]  •  0  THEN  U.DF[I]:-0  ELSE  U.DF[I]:=-G.DF[I]*U.F/G.F; 
RES:*U 

END; 

OPERATOR  /  (G:  GRADIENT ;R:  REAL)  RES:  GRADIENT; 

VAR  U:  GRADIENT; I:  OIMTYPE; 

BEGIN  U.F:*G.F/R;FOR  I:-1  TO  DIM  DO  U.DF[I]:-G.DF[I]/R;RES:-U  END; 


'  '•* 


# 


OPERATOR  /  (GA:  GRADIENT ;GB:  GRADIENT)  RES:  GRADIENT; 

VAR  U:  GRADIENT ;I :  DIKTYPE; 

BEGIN  U.F:«GA.F/GB.F;FOR  I:»l  TO  DIM  DO 

BEGIN  IF  GA.DF[I]  -  0  THEN  U.DF[I]:«0  ELSE  U.DF[I]:=GA.DF[I]/GB.F; 
IF  GB.DF[I]  <>0  THENU.DF[I]:*U.DF[I]-GB.DF[I]*U.F/GB.F 

END; 

RES:-U 

END; 


Power  Operators  (File  GRAD_POW . S ) 

OPERATOR  **  (R:  REAL ;K :  INTEGER)  RES:  REAL; 

VAR  L:  INTEGERS:  REAL; 

BEGIN  IF  K  < ■  0  THEN  Ur-1/R; 

IF  K  =  0  THEN  U:=l; 

ELSE  IF  K  -  1  THEN  U:=R 

ELSE  BEGIN  L:»ABS(K) ;U:«1 ;REPEAT  IF  L  MOD  2  *  1 
THEN  U:«R*U;L:«L  DIV  2;  IF  L  O0 
THEN  R:-R*R  UNTIL  L  -  0; 

IF  K  (0  THEN  U:*l/U 

END; 

RES:=U 

END; 

OPERATOR  **  (RA:  REAL;RB:  REAL)  RES:  REAL; 

VAR  K:  INTEGER ;U:  REAL; 

BEGIN  IF  RA  =  1  THEN  U:*l  ELSE  IF  (RA-0)  AND  (RB  > 0)  THEN  U:*0  ELSE 
BEGIN  K:»TRUNC(RB);U:«RA**K;IF  RB  OK  THEN 
U:=U*EXP{ (RB-K)*LN(RA) ) ; 


ENO; 


OPERATOR  **  (G:  GRADIENT ;K:  INTEGER)  RES:  GRADIENT; 

VAR  I:  DIMTYPE;M:  REAL;U:  GRADIENT; 

BEGIN  U.F:=G.F**K;IF  K  *  1  THEN  U.DF:»G.DF  ELSE 

IF  K  *  0  THEN  FOR  I:=l  TO  DIM  DO  U.DF[I]:«0  ELSE 
IF  (G.F-0)  AND  {K  >1)  THEN  FOR  I:-l  TO  DIM  DO  U.DF[1]:*0  ELSE 
BEGIN  M:*K*U.F/G.F;FOR  I:-l  TO  DIM  DO  IF  G.DF[I]  -  0  THEN 
U.DF[I]:=0  ELSE  U.DF[I]:-M*G.DF[I] 

END; 

RES.-U 

END; 

OPERATOR  **  (G:  GRADIENT;R:  REAL)  RES:  GRADIENT; 

VAR  I:  DIMTVPE;M:  REAL;U:  GRADIENT; 

BEGIN  U.F:=G.F**R;IF  R  -  I  THEN  U.DF:=G.DF  ELSE 

IF  R  -  0  THEN  FOR  I:«1  TO  DIM  DO  U.DF[I]:-0  ELSE 
IF  (G.F-0)  AND  (R  >1)  THEN  FOR  I:-l  TO  DIM  DO  U.DF[I]:»0  ELSE 
BEGIN  M:-R*U.F/G.F;FOR  I:-l  TO  DIM  DO  IF  G.DF[I]  *  0  THEN 
U.DF[I]:«0  ELSE  U.DF[I] :»M*G.DF[I] 

END; 

RES:»U 

END; 

OPERATOR  **  (R:  REAL;G:  GRADIENT)  RES:  GRADIENT; 

VAR  I:  DIMTYPE;M:  REAL;U:  GRADIENT; 

BEGIN  U.F:»R**G.F;IF  (R»0)  AND  (G.F  >0)  THEN  FOR  I:»l  TO  OIM  00 
U.DF[I]:*0  ELSE 

BEGIN  M:-U.F*LN(R);FOR  I:»1  TO  DIM  DO  IF  G.DF[I]  -  0  THEN 
U.DF[I]:*0  ELSE  U.DF[I]:-M*G.DF[I] 

END; 

RES:-U 

END; 

OPERATOR  **  (K:  INTEGER;G:  GRADIENT)  RES:  GRADIENT; 

VAR  R:  REAL;U:  GRADIENT; 

BEGIN  R:-K;U:-R**G; 


OPERATOR  **  (GA:  GRADIENT ;GB:  GRADIENT)  RES:  GRADIENT; 
VAR  I:  DIMTYPE;R:  REAL;U,V:  GRADIENT; 

BEGIN  R :  -GA .  F  ;U : «R**G8 ;R : "GB . F ; V : »GA**R ; 

FOR  I :■!  TO  OIM  DO  U.DF[I].-U.DF[l]+V.DF[I] 
RES:*U 

END; 
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APPENDIX  B.  SOURCE  CODE  FOR  STANDARD  GRADIENT  FUNCTIONS 


(File  GRAD  FUN  .S) 


B.l.  Absolute  Value 

FUNCTION  GABS(G:  GRADIENT):  GRADIENT; 

VAR  I:  DIMTYPE ;M:  REAL ;U:  GRADIENT; 

BEGIN  U.F:*ABS(G.F);M:*G.F/U.F;FOR  I :-l  TO  DIM  DO  IF  G.DF[I]  -  0  THEN 
U.DF[I]:«0  ELSE  U.DF[I]:-M*6.DF[I]; 

GABS:*U 

END; 

B.2.  Square  Root 

FUNCTION  GSQRT(6:  GRADIENT):  GRADIENT; 

VAR  I:  DIMTYPE ;M:  REAL ;U :  GRADIENT; 

BEGIN  U.F:-SQRT(G.F);M:-2*U.F;F0R  I:-l  TO  DIM  DO  IF  G.DF[I]  -  0 
THEN  U.DF[I]:-0  ELSE  U.DF[I]:»G.DFf I]/M; 

GSQRT :*U 

END; 

B.3.  Exponential  (Base  e) 

FUNCTION  GEXP(G:  GRADIENT):  GRADIENT; 

VAR  I:  DIMTYPE ;U:  GRADIENT; 

BEGIN  U.F:-EXP(G.F);FOR  I:«1  TO  DIM  DO  IF  G.DF[I]  «  0  THEN  U.DF[I]:-0 
ELSE  U.DF[I]:-U.F*G.DF[I]; 

GEXP:*U 

END; 

B.4.  Natural  Logarltha 

FUNCTION  GLN(G:  GRADIENT):  GRADIENT; 

VAR  I:  DIMTYPE ;U:  GRADIENT; 

BEGIN  U.F:»LN(G.F);FOR  I:-l  TO  DIM  DO  IF  G.DF[I]  -  0  THEN  U.DF[I]:«0 
ELSE  U.OF[I]:«G.DF[I]/G.F; 

GLN:*U 
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B.5.  Arctangent 


FUNCTION  GARCTAN(G:  GRADIENT):  GRADIENT; 

VAR  I:  DIHTYPE;M:  REAL;U:  GRADIENT; 

BEGIN  U.F:«ARCTAN(G.F);M:«l+G.F*G.F;FOR  I:-l  TO  DIH  DO  IF  G.DF[I] 
THEN  U.DF[I]:-0  ELSE  U.DF[I]:-G.DF[I]/M; 

GARCTAN:=U 

END; 


B.6.  Sine 


FUNCTION  GSIN(G:  GRADIENT):  GRADIENT; 

VAR  I:  DIMTYPE;N:  REAL;U:  GRADIENT; 

BEGIN  U.F:«SIN(G.F);M:«COS(G.F);FOR  I:-l  TO  DIN  DO  IF  G.DF[I]  -  9 
THEN  U.DF[I]:*0  ELSE  U.DF[I]:-H*G.DF[I]; 

GSIN:*U 


FUNCTION  GCOS(G:  GRADIENT):  GRADIENT; 

VAR  I:  DINTYPE;M:  REAL;U:  GRADIENT; 

BEGIN  U.F:*COS(G.F);N:»-SIN(G.F);FOR  I:»l  TO  DIN  DO  IF  G.DF[I]  «  9 
THEN  U.DF[I]:*0  ELSE  U.DF[I]:=H*G.DF[I] 

END; 


APPENDIX  C.  OUTPUT  OF  THE  GRADIENT  PROGRAM  FOR  TOE 
SOLUTION  OF  SYSTEM  (9.10)  BY  NEWTON'S  METHOD 


INITIAL  VALUES 

X  ■  I.OOOOOOOOOOOE+OO 
X  ■  1 .OOOOOOOOOOOE+OO 
Z  *  1.00000000000E+00 
FUNCTION  VALUES 

F(X,V,Z)  -  1 . 70000000000E+01 
fi(X.V,Z)  ■  O.OOOOOOOOOOOE+OO 
H(X,V,Z)  -  O.OOOOOOOOOOOE+OO 
CORRECTION  VECTOR 

OX  «  -7 . 08333333334E -02 
DV  »  -2.1 2500000000E -01 
OZ  «  2 . 83333333334E -01 

VALUES  FROM  ITERATION  NUMBER  1 
X  -  9.29166666667E-01 
Y  •  7. 8750000 00 OOE-Ol 
Z  »  1 .28333333333E+00 
FUNCTION  VALUES 

F(X.V.Z)  -  4 . 791 91 71 3930E+00 
G(X,Y,Z)  =  1. 30451 388890E-01 
H(X,Y,Z)  »  1 .4696686 9220E-02 
CORRECTION  VECTOR 

OX  -  -4.20921 371 897E -02 
DY  -  -9.43241406975E-02 
OZ  -  3.7531 3068779E -02 

VALUES  FROM  ITERATION  NUMBER  2 
X  «  8.87074529477E-01 
Y  -  6.93I75859303E-01 
Z  ■  1. 3208646402 1E+00 


ITERATION  NUMBER  2  (CONTINUED) 
FUNCTION  VALUES 

F(X,Y,Z)  -  6 .45309522200E-01 
<»(X,Y,Z)  *  1 .20773905300E-02 
H(X,Y,Z)  *  4 . 8641 7092500E-03 
CORRECTION  VECTOR 

DX  *  -8.8301 31 1 3461 E-03 
OY  -  -1.59811519852E-02 
DZ  =  9.74516049898E-03 

VALUES  FROM  ITERATION  NUMBER  3 
X  *  8.78244398342E-01 
Y  «  6.77194707318E-01 
Z  -  1. 33060980071 E+00 
FUNCTION  VALUES 

F(X,Y.Z)  -  1 .84509451000E-02 
G(X,Y,Z)  -  4 . 28336590000E-04 
H(X,Y,Z)  -  2 . 0681 034 1 OOOE-04 
CORRECTION  VECTOR 

OX  -  -2.78405682738E-04 
DY  -  -4. 3740361 2547E-04 
DZ  -  2.4541 1801 068E-04 

VALUES  FROM  ITERATION  NUMBER  4 
X  -  8 . 77965992659E-01 
Y  ■  6 . 76757303705E-01 
Z  *  1.33085521 251 E+00 
FUNCTION  VALUES 

F(X,Y,Z)  *  1 .47972000000E-05 
G(X,Y,Z)  *  3. 290500 0000OE-07 
H(X,Y,Z)  «  2 . 041 96000000E-07 


ITERATION  NUMBER  4  (CONTINUED) 

CORRECTION  VECTOR 

DX  »  -2. 3238361 7048E -07 
DY  *  -3 . 331 84805859E -07 
DZ  ■  1.99108934899E-07 

VALUES  FROM  ITERATION  NUMBER  5 
X  -  8.77965760275E-01 

Y  *  6.76756970520E-01 
Z  =  1. 33085541 162E+00 

FUNCTION  VALUES 

F(X,Y,Z)  *  1.00000000000E-10 
G(X,Y,Z)  •  O.OOOOOOOOOOOE+OO 
H(X,Y,Z)  * -1.00000000000E-12 

CORRECTION  VECTOR 

DX  «  -1 .181 97463333E-12 
DY  =  -3 . 73328280534E-1 2 
DZ  «  2.6781 71 03788E-1 2 

VALUES  FROM  ITERATION  NUMBER  6 
X  =  8.77965760274E-01 

Y  *  6.76756970S16E-01 
Z  *  1. 33085541 162E+00 

FUNCTION  VALUES 

F{X,Y,Z)  *  O.OOOOOOOOOOOE+OO 
G(X,Y,Z)  *  O.OOOOOOOOOOOE+OO 
H(X,Y,Z)  =  2 . OOOOOOOOOOOE - 1 2 

CORRECTION  VECTOR 

DX  «  -4.1 855801 9594E-1 3 
DY  ■  1.03209645475E-12 

DZ  »  -2.4871 1360538E-1 3 


VALUES  FROM  ITERATION  NUMBER  7 
X  *  8.77965760274E-01 

Y  *  6.7675697051 7E-0I 
Z  -  1. 330855411 62E+00 

FUNCTION  VALUES 

F(X,Y,Z)  ■  O.OOOOOOOOOOOE+OO 
G(X,Y,Z)  *  O.OOOOOOOOOOOE+OO 
H(X,Y,Z)  *  1. OOOOOOOOOOOE- 12 

CORRECTION  VECTOR 

DX  =  -2 . 09279009797E- 1 3 
DY  *  5 .16048227375E-1 3 
DZ  «  -1.24355680269E-13 

VALUES  FROM  ITERATION  NUMBER  8 
X  *  8.77965760274E-01 

Y  =  6. 7675697051 8E -01 
Z  *  1. 3308554 1162E+00 

FUNCTION  VALUES 

F(X,Y,Z)  =  O.OOOOOOOOOOOE+OO 
G(X,Y,Z)  =  O.OOOOOOOOOOOE+OO 
H(X,Y,Z)  =  O.OOOOOOOOOOOE+OO 

CORRECTION  VECTOR 

DX  =  O.OOOOOOOOOOOE+OO 
DY  =  O.OOOOOOOOOOOE+OO 
DZ  -  O.OOOOOOOOOOOE+OO 
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